home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / multiroots / newton.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-12-09  |  3.7 KB  |  152 lines
  1. /* multiroots/newton.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21.  
  22. #include <stddef.h>
  23. #include <stdlib.h>
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <float.h>
  27.  
  28. #include <gsl/gsl_math.h>
  29. #include <gsl/gsl_errno.h>
  30. #include <gsl/gsl_multiroots.h>
  31. #include <gsl/gsl_linalg.h>
  32.  
  33. typedef struct
  34.   {
  35.     gsl_matrix * lu;
  36.     gsl_permutation * permutation;
  37.   }
  38. newton_state_t;
  39.  
  40. static int newton_alloc (void * vstate, size_t n);
  41. static int newton_set (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  42. static int newton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  43. static void newton_free (void * vstate);
  44.  
  45. static int
  46. newton_alloc (void * vstate, size_t n)
  47. {
  48.   newton_state_t * state = (newton_state_t *) vstate;
  49.   gsl_permutation * p;
  50.   gsl_matrix * m;
  51.  
  52.   m = gsl_matrix_calloc (n,n);
  53.   
  54.   if (m == 0) 
  55.     {
  56.       GSL_ERROR_VAL ("failed to allocate space for lu", GSL_ENOMEM, 0);
  57.     }
  58.  
  59.   state->lu = m ;
  60.  
  61.   p = gsl_permutation_calloc (n);
  62.  
  63.   if (p == 0)
  64.     {
  65.       gsl_matrix_free(m);
  66.  
  67.       GSL_ERROR_VAL ("failed to allocate space for permutation", GSL_ENOMEM, 0);
  68.     }
  69.  
  70.   state->permutation = p ;
  71.  
  72.   return GSL_SUCCESS;
  73. }
  74.  
  75. static int 
  76. newton_set (void * vstate, gsl_multiroot_function_fdf * FDF, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
  77. {
  78.   newton_state_t * state = (newton_state_t *) vstate;
  79.  
  80.   size_t i, n = FDF->n ;
  81.  
  82.   state = 0 ; /* avoid warnings about unused parameters */
  83.  
  84.   GSL_MULTIROOT_FN_EVAL_F_DF (FDF, x, f, J);
  85.  
  86.   for (i = 0; i < n; i++)
  87.     {
  88.       gsl_vector_set (dx, i, 0.0);
  89.     }
  90.  
  91.   return GSL_SUCCESS;
  92. }
  93.  
  94. static int
  95. newton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
  96. {
  97.   newton_state_t * state = (newton_state_t *) vstate;
  98.   
  99.   int signum;
  100.  
  101.   size_t i;
  102.  
  103.   size_t n = fdf->n ;
  104.  
  105.   gsl_matrix_memcpy (state->lu, J);
  106.  
  107.   gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
  108.  
  109.   gsl_linalg_LU_solve (state->lu, state->permutation, f, dx);
  110.       
  111.   for (i = 0; i < n; i++)
  112.     {
  113.       double e = gsl_vector_get (dx, i);
  114.       double y = gsl_vector_get (x, i);
  115.       gsl_vector_set (dx, i, -e);
  116.       gsl_vector_set (x, i, y - e);
  117.     }
  118.  
  119.   {
  120.     int status = GSL_MULTIROOT_FN_EVAL_F_DF (fdf, x, f, J);
  121.     
  122.     if (status != GSL_SUCCESS) 
  123.       {
  124.         return GSL_EBADFUNC;
  125.       }
  126.   }
  127.  
  128.   return GSL_SUCCESS;
  129. }
  130.  
  131.  
  132. static void
  133. newton_free (void * vstate)
  134. {
  135.   newton_state_t * state = (newton_state_t *) vstate;
  136.  
  137.   gsl_matrix_free(state->lu);
  138.  
  139.   gsl_permutation_free(state->permutation);
  140. }
  141.  
  142.  
  143. static const gsl_multiroot_fdfsolver_type newton_type =
  144. {"newton",                /* name */
  145.  sizeof (newton_state_t),
  146.  &newton_alloc,
  147.  &newton_set,
  148.  &newton_iterate,
  149.  &newton_free};
  150.  
  151. const gsl_multiroot_fdfsolver_type  * gsl_multiroot_fdfsolver_newton = &newton_type;
  152.